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We study the unitary time evolution of antiferromagnetic order in the Hubbard model after a 
quench starting from the perfect Neel state. In this setup, which is well suited for experiments 
with cold atoms, one can distinguish fundamentally different pathways for melting of long-range 
order at weak and strong interaction. In the Mott insulating regime, melting of long-range order 
occurs due to the ultra-fast transfer of energy from charge excitations to the spin background, while 
local magnetic moments and their exchange coupling persist during the process. The latter can be 
demonstrated by a local spin-precession experiment. At weak interaction, local moments decay along 
with the long-range order. The dynamics is governed by residual quasiparticles, which are reflected 
in oscillations of the off-diagonal components of the momentum distribution. Such oscillations 
provide an alternative route to study the prethermalization phenomenon and its influence on the 
dynamics away from the integrable (noninteracting) limit. The Hubbard model is solved within 
nonequilibrium dynamical mean-field theory, using the density matrix-renormalization group as an 
impurity solver. 

PACS numbers: 71.10.Fd,75.40.Mg 


I. INTRODUCTION 

Ultra-fast pump-probe experiments on condensed mat¬ 
ter systems and experiments with cold gases in optical 
lattices have opened the intriguing possibility of control¬ 
ling transitions between complex phases on microscopic 
timescales. This has motivated intensive theoretical ef¬ 
forts to understand fundamental aspects of the dynamics 
in interacting many-body systems, and leads to predic¬ 
tions in marked contrast to the naive expectation that 
interactions imply rapid thermalization [l|: Integrable 
systems can keep memory of the initial state for all times 
and relax to a generalized Gibbs ensemble ii , but also 
away from integrabiliW thermalization can be delayed by 
prethermalizaton |j-l7|, and one can identify regimes of 
different dynamical behavior which are clearly separated 
by non-thermal critical points dSl. 

Of particular interest with respect to complex phases 
in condensed matter is the dynamics of symmetry broken- 
states dm. While the relevant relaxation mechanisms 
after a perturbation are hard to disentangle in a solid, 
cold atoms in optical lattices provide a versatile platform 
to investigate isolated quantum systems in ideal situa¬ 
tions. The preparation of thermodynamic long-range or¬ 
dered phases in cold atoms is still a challenge d3 : but ad¬ 
vanced techniques for lattice design have made it possible 
to prepare an ordered state on a lattice of isolated sites. 
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and to probe its dynamics after tunneling between the 
sites is switched on [TBUlj . In the following we consider 
such a setup for the Fermi-Hubbard model, a paradigm 
model for emergent long-range order in condensed mat¬ 
ter systems. We will simulate the time evolution starting 
from a classical Neel state in which neighboring lattice 
sites of a bipartite lattice are occupied with particles of 
opposite spin. 

In general, one can anticipate fundamentally differ¬ 
ent pathways for melting of long-range antiferromagnetic 
order in the weakly and strongly interacting Hubbard 
model: For strong interaction, long-range order arises 
from antiferromagnetically coupled local moments, which 
emerge when charge fluctuations are frozen. Magnetic 
order could thus possibly melt via the destruction of the 
local moments themselves, through a reduction of the ef¬ 
fective exchange interaction (while moments persist), 
or, along a quasi-thermal pathway, by the transfer of en¬ 
ergy from excited quasiparticles (hot electrons) to spins. 
The latter mechanism is intensively studied in the con¬ 
text of photo-carrier relaxation in high-Tc cuprates (2^ 
[ 2 ^, where the investigation of the spin-charge interac¬ 
tion challenges the limits for the time-resolution in state 
of the art pump-probe experiments (29l - [^ . For weak 
interaction, on the other hand, quasiparticle states may 
be important to understand relaxation processes. In the 
paramagnetic phase the conservation of the quasiparticle 
momentum occupations imposes constraints on the dy¬ 
namics which can lead to prethermalization sai& 
[ 3 ^. Prethermalizaton, which was recently observed in 
a one-dimensional Bose gas Q, has been suggested to 
be a universal feature of near-integrable systems Q , but 
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previous predictions for the Hubbard model rely on a dis¬ 
continuity of the momentum distribution which is absent 
at nonzero temperature and thus experimentally hard to 
observe. Here we show that the symmetry-broken initial 
state provides an alternative perspective to investigate 
this physics and its breakdown far from integrability. 

Quenches from a Neel state have been explored in 
quantum spin models also as a way to prepare 

ordered states in the Hubbard model [s^ , but a pure spin 
model cannot describe the relevant dynamics of charge 
excitations and local moments. The Hubbard model has 
been studied in one dimension u sing the density-matrix 
renormalization group (DMRG) [39|. For the dynamics 
of lattice fermion models in more than one dimension, 
nonequilibrium dynamical mean-field theory (DMFT) 
is the most promising approach. Quenches within 
the antiferroma gne tic phase of the Hubbard model at 
strong-coupling |24l | are in line with the “quasi-thermal” 
pathway discussed above. The regime of intermediate in¬ 
teractions, where the notion of local moments becomes 
ambiguous, or to weak coupling, where prethermaliza- 
tion may be expected, has been elusive so far. Previous 
numerical solutions of the DMFT equations were based 
on the self-consistent strong-coupling expansion [4l| or 
weak-coupling impurity solvers |13l l42l . l43l | , which both 
fail at intermediate coupling, while weak-coupling quan¬ 
tum Monte Carlo studies [a, |4l| are most efficient for 
noninteracting initial states and restricted to short times. 
In this work we overcome these limitations using a re¬ 
cently developed Hamiltonian based formulation for the 
impurity model of nonequilibrium DMFT (3 | , which has 
opened the possibility to use wave-function based tech¬ 
niques to solve the DMFT equations [1^ Here we 
use DMRG as an impurity solver , which allows us to 
reach sufficiently long times in the evolution to address 
the above issues. 


II. MODEL AND METHODS 

Throughout this work we consider the single-band 
Hubbard model at half-filling, with nearest-neighbor hop¬ 
ping J and on-site Coulomb repulsion U. The Hamilto¬ 
nian is given by 

H = -J{t) ^ + (1) 

where (cia-) are electron creation (annihilation) op¬ 
erators for lattice site i and spin a, and riia = cj^Cicr. 
The model is solved using nonequilibrium DMFT [40l |. 
for a Bethe lattice in the limit of infinite coordination 
number Z and hopping J = ^ where the approach 

becomes exact [43 . The energy unit is set by J* = 1, and 
time is measured in inverse energy, i.e, the free density of 
states is given by D{e) = \/4 — e2/(27r). To simulate the 
quench, we choose a time-dependent hopping J*(t) = 0 
for t < 0 and J*(t) = 1 for t > 0. For t <0, the system 


therefore consists of a set of isolated lattice sites, which 
are prepared in a classical Neel state, 

|'kNeel) = n4tnS4|0)^ (2) 

ieA jgB 

where A and B are sub-lattices of the bipartite Bethe 
lattice. 

In DMFT, the lattice model is mapped to a set of im¬ 
purity problems, one for each inequivalent lattice site 
j = A,B, with a time-dependent hybridization func¬ 
tion Ajcr(t) t'). (In this expression, time arguments lie on 
the Keldysh contour; see Ref. [i^ for a detailed descrip¬ 
tion of nonequilibrium DMFT and the Keldysh formal¬ 
ism.) For the Bethe lattice, the latter is determined self- 
consistently by 

where = ~*(TcCj(T(t)c|^(t')) is the local 

Green’s function. To solve the impurity model with a non 
time-translationally invariant hybridization function, we 
derive an equivalent representation in terms of a time- 
dependent Anderson impurity Hamiltonian [i^l with up 
to L = 24 bath orbitals, from which the time-dependent 
Green’s functions are computed using a Krylov time- 
propagation for matrix product states [d^. The Hamil¬ 
tonian representation of the DMFT impurity model is 
exact for small times, but an increasii^ number of bath 
sites is needed to reach longer times [^. We verify the 
convergence of the solution with the bath size L. Up to 
L = 12, the results have also been cross checked with 
a Krylov time-propagation in the full Hilbert space. For 
further details of the numerical solution, see Appendix lAl 


III. RESULTS 

Figure [T] shows the time evolution of the antiferromag¬ 
netic order parameter M{t) and the double occupation 
d{t) = {n^{t)ni{t)) after the quench, for various values of 
the Coulomb interaction. In order to account for the triv¬ 
ial reduction of the local spin expectation value by vir¬ 
tual charge fluctuations, we define M{t) as the staggered 
order M^tagg = (riAtit) - riAiit)) = (riBiit) - UBtit)), 
normalized by the probability Pi for a site to be singly 
occupied, Pi{t) = 1 — 2d{t). To test whether the system 
thermalizes after the quench, we compare to an equilib¬ 
rium state at the same internal energy (which is zero for 
the Neel state). The corresponding effective temperature 
Teff (Fig.[TJ:) lies above the Neel temperature TNeei for all 
values of U . This implies a paramagnetic state after 
thermalization. While M{t) indeed continues to decay 
throughout the simulated time interval, the double occu¬ 
pancy saturates to a non-thermal value for C/ > 4 (ar¬ 
rows in Fig. [T}d point at the thermalized value d(Teff)), 
in agreement with earlier studies on the lifetime of dou- 
blons in the paramagnetic Mott regime [sTI - fs^ . Already 
at a first glance, the relaxation of M{t) and d{t) there¬ 
fore suggests different mechanisms for small and large 
values of 17, with a rapid and oscillatory decay of M{t), 
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FIG. 1: a) Time evolution of the order parameter M{t) = 
A^stagg/[1 —2d(t)] for different values of the Coulomb repulsion 
in the range 17 = 0 to 17 = 10. b) The double occupation d{t) 
at the same values of U. Arrows indicate the double occu¬ 
pation in a thermalized state at the same total energy as the 
quenched state (as obtained from equilibrium DMFT, using 
Continuous-time Quantum Monte Carlo as an impurity 
solver). All thermalized states are in the paramagnetic phase, 
the corresponding inverse temperature 1 /Teff is plotted in c). 


and a long-lived non-thermal state, respectively. In the 
following we will analyze the two regimes in more detail. 


A. Weak-coupling: Residual quasiparticles 

For quenches to small U the Hamiltonian is close to the 
integrable point U = Q. This suggests to study relaxation 
in terms of the momentum occupation nk(t) = {c\ck), 
which is conserved at t/ = 0. For a state with trans¬ 
lational symmetry breaking, the single-particle density 
matrix pkk'{t) = (cjc/(t)cfc(t)) is no longer diagonal in 
momentum k. (The discussion holds for a general lattice 
like the Bethe lattice when k denotes the eigenstates of 
the translationally invariant hopping matrix.) For near¬ 
est neighbor hopping on a bipartite lattice, eigenstates 
come in pairs k,k with single particle energy = — e^, 
where the wave functions for k and k differ by a stag¬ 
gered phase = ± for f G A{B), and Pkk 7 ^ 0 if the 
symmetry between sub-lattices is broken. (On the cu¬ 
bic lattice, k and k = k + [it, tt, ...) are momenta re¬ 
lated by the antiferromagnetic nesting vector.) In Fig. [2] 
we plot the diagonal and off-diagonal components of the 
single-particle density matrix in terms of the two func¬ 
tions n{ek,t) = {cl{t)ck{t)) and TT{ek,t) = {c\{t)ck{t)), 
which depend on k only via due to the locality of 
the self-energy within DMFT. In the thermalized state. 


7r(e) = 0 because the state does not break the sub¬ 
lattice symmetry, while in the localized initial state, 
n{e,t = 0) = 7r(e,t = 0) = 1/2. In agreement with the 
behavior of the double occupancy, n(e, t) does not ther- 
malize at large U (thermalized values and 

are shown by solid lines). For U < 2, however, differ¬ 
ences between n(e,t) and nT^jj{e) become tiny. This is in 
stark contrast to the behavior of the paramagnetic sys¬ 
tem after a quench from [7 = 0, where prethermalization 
manifests itself precisely in the difference between n(e, t) 
and nTgif{e) Moreover, around U = 3 the relaxation 
of 7r(e,t) changes from an oscillatory to a monotonous 
decay. (In Fig. [2 we plot the real part of 7r(e,t), the 
imaginary part shows a similar crossover from oscillatory 
to non-oscillatory behavior.) 

This observation may be explained following the per¬ 
turbative arguments of Refs, jl, @|. To second order in 
U, the Hamiltonian m is unitarily equivalent to a model 
H = J2k ^kCkCk + 0(U^) which is quadratic in terms of 
quasiparticle operators Ck', we have Ck = RkCk + incoh. 
with a finite residue Rk, where incoh. denotes incoher¬ 
ent contributions, i.e., an admixture of particle-hole ex¬ 
citations to higher order in U. Hence the momentum 
occupation is given by nk{t) = R^{c\.{t)ck{t)) -I- incoh.. 
The term proportional to R\ (the coherent part) is un¬ 
changed by the time evolution to second order in U. 
Back-transforming gives rifc(t) = Rf,nk{0)+incoh., where 
the incoherent contribution is a smooth function of k. For 
quenches in the paramagnetic phase, nk{t) thus preserves 
the initial discontinuity at the Fermi surface, which can 
be taken as a measure of prethermalization Q. In the 
symmetry broken state, however, nk(0) is independent 
of k, and thus nk(t) does not clearly exhibit the exis¬ 
tence of residual quasiparticles. In fact, the numerical 
results suggest that the incoherent part can accurately 
be described by a thermal distribution. In contrast, a 
similar argument for the off-diagonal component shows 
that 7rk(t) = Rk(cl(t)ck(t)} + incoh. = t = 

0 ) + incoh., where we have used the time evolution of 
the quasiparticle, Ck(k){t) = e^*'^''*C;.(^)(0) and Rk = Rk- 
Hence we find that the residual quasiparticle dynamics 
leading to prethermalization close to the integrable point 
U = 0 can be studied very conveniently with the sym¬ 
metry broken initial state in terms of oscillations in the 
off-diagonal components of the momentum occupation. 

Similar to the interaction quench in the paramagnetic 
phase iini , we find that the “prethermalization” regime 
in which residual quasiparticles dominate the dynamics is 
limited to small interactions; at large interactions, 7r(e, t) 
relaxes to zero monotonously (see the U = 6 data in 
Fig. [2b) and the distribution becomes flat over the Bril- 
louin zone. Below we will see that the dynamics at large 
U can be analyzed in terms of well defined localized mo¬ 
ments. In contrast to the quench in the paramagnetic 
phase, the crossover between the weak and strong cou¬ 
pling regimes is relatively smooth and occurs between 
U = 2 and U = 3: In Fig. |3b, we exemplarily plot 
Re7r(e,t) for fixed e = —1.5 and various U. For 17 < 2, 
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FIG. 2: Diagonal component n{ek,t) = (c|.Cfc) of the momentum occupation (panels a, c, and e) and off-diagonal component 
Re7r(efc) = (cj.cj.) (panels b, d, and f), where k and k are pairs of single particle states coupled by a staggered potential, plotted 
for quenches to three different values of 17 as a function of the energy eu ranging from —2 to 2 in the band of the Bethe lattice. 
The bold lines indicate momentum distributions obtained in the (paramagnetic) equilibrium state at the same energy. 
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FIG. 3: a) Line-out of for t = —1.5 and vari¬ 

ous values of U. Solid lines are fits with the sum of a 
decaying exponential background and decaying oscillations, 
f{t) — 6exp(—c(t — to)) + aexp(—r(t — to)) cos(2e't -|- (ji) for 
t > to = 1.5. b) Amplitudes of the background, b, the oscil¬ 
lations, a, the quasi-particle energy e and the quasi-particle 
decay rate F as a function of U. 


the curves can be accurately fit with decaying oscilla¬ 
tions /i(t) = aexp(—Ft) cos(—2e't -I- (/>), where in agree¬ 
ment with the discussion above the quasi-particle en¬ 
ergy e' —>■ e, and F ^ 17^ for [/ —>■ 0 (solid lines in 


Fig. [3^, fit parameters in Fig. [3}d). For 17 > 3, on the 
other hand, a good fit is a monotonously decaying curve 
/ 2 (t) = l)exp(—ct). For 2 < 17 < 3, there is a crossover 
between the two behaviors, as evidenced by the depen¬ 
dence of the amplitudes a and b of the monotonous and 
the oscillating component on U (Fig. [5 )d). 

Before discussing the strong-coupling regime, we note 
that off-diagonal momentum distributions can in princi¬ 
ple be measured by a modified time-of-flight measure¬ 
ment, if before releasing the cloud, one would switch 
off the tunneling and the interaction, switch on a stag¬ 
gered potential which is ±A on the A and B sub-lattice, 
respectively, and evolve for a given time tm- Time of 
flight measures the regular momentum occupation Uk = 
{c\ck) after that procedure. In the k, k basis, the stag¬ 
gered potential is given by 77 = A^^(c^c^ -I- H.c.), so 
that nk{t + tm) = nfe(t) cos^(tmA)nfe(t) -I- sin^(tmA) -|- 
sin(2tmA)Im7rfe(l:) after propagation in the pure stag¬ 
gered potential from time t to t -\- trm and Im7rfe(t) can 
be extracted. 


B. Dynamics of local moments 

In a Mott insulator at large U one can expect the ex¬ 
istence of well-defined local moments. It is an intriguing 
question whether these moments persist in the quenched 
state while the long-range order disappears, and to what 
extent the crossover in relaxation behavior from weak to 
strong coupling can be characterized in terms of these 
local moments. In the following, we propose a simple 
experiment to distinguish the existence and strength of 
moments in the quenched state: one spin in the initial 
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FIG. 4: Dynamics at a probe site o, where the spin is initially flipped to the x-direction. a)-c) Expectation values {Sz{t)), 
{Sy{t)), and {Sx{t)) for various values of U (see legend in d)). d) Trajectory of the spin (S) in the Sx-Sy plane. The inset 
illustrates the initial state, with one spin flipped to the x direction, and the effective exchange field, d) Effective exchange 
interaction, d(()(f)/df/|Mstagg|, where (f>{t) = atan(5'j;/5'a,) is the angle in the Sx-Sy plane, for U = 3,4,5,7,10. The dotted 
lines correspond to the perturbative value of the exchange interaction Jex = 4 J* /U. f) Temperature-dependent local moment 
in equilibrium, defined by (1//I) dT(5'z(r)S'z(0)), obtained using CTQMC as an impurity solver. 


Neel state on a given site (the probe site “o”) is flipped 
to the x-direction (see Fig. |4ji, inset). Choosing o on the 
A-sublattice of the Neel state, the initial state ([2]) of the 
dynamics is changed to(|'FNeei,t) + l^ Neel) 4-) ) / , where 
jllfNeebO-) = cj ^Co't-l^'Neei)- In & perfect local moment 
picture, the spin should then precess in the exchange field 
of it’s neighbors. 

The inhomogeneous setup with one probe spin can be 
solved within DMFT, where it corresponds to a modified 
impurity problem at site o, while the rest of the lattice 
is unchanged (see App. |^. Figure H^-c shows the local 
spin expectation values (Sz), (Sy), and (Sx) at site o for 
various values of the interaction. In Fig. we show the 
trajectory of the spin in the Sx-Sy plane, starting from 
S'a; = 1, S';; = 0 at time t = 0. For large U one can indeed 
observe a precessional motion in the Sx-Sy plane, as ex¬ 
pected for a local moment subject to an exchange field 
in the z-direction. For [/ = 0, on the other hand, the 
spin-dynamics is entirely longitudinal, showing no sign 
of well-defined local moments. (For {7 = 0, the dynam¬ 
ics can be solved analytically, yielding Sx,o = 
while So,y = 0 for the Bethe lattice at Z = oo, where 
Ji{x) is the first Bessel function (Appendix [C])). There is 
a crossover between the two relaxation regimes. 

Although the exchange interaction is in principle not 
an instantaneous interaction on the timescale of the elec¬ 
tronic hopping (2^, it is illustrative to quantify the 
precession dynamics in terms of an effective exchange 
field. For this purpose we follow Refs. [H, and de¬ 
fine BeS such that {S{t)) satisfies the equation of mo¬ 
tion d/dt{S{t)) = Beff X {S{t)). We can assume that 
BeS = BeffZ acts only in the z direction (parallel to the 


order parameter Mgtagg on the neighboring sites), and 
use the parametrization B^s = dexAfstagg to define an ef¬ 
fective exchange interaction Jex; the latter is then given 
by Jex = (^(t)/|Mstagg| where 4>{t) = ataa{Sy{t) / Sx{t)) 
is the angle of the spin in the x-y plane. The resulting 
value Jex is plotted in Fig. |3^. For large U, Jex shows 
very good agreement with the perturbative value of the 
exchange in the Hubbard model, 4J^/{7, and is not sub¬ 
stantially decreasing with time even for quenches at inter¬ 
mediate interaction ({7 ~ 4) where the order parameter 
quickly decays to zero (see Fig. [T]). Equilibrium esti¬ 
mates of the local moment in the intermediate coupling 
regime (Fig.[3f) furthermore show tendencies of moment 
formation at elevated temperatures, which may explain 
why some spin precession occurs even for {7 = 2. The 
combination of these results shows that the melting of 
long-range order proceeds by the “quasi-thermal” path¬ 
way discussed in the introduction, i.e., a disordering of 
exchange-coupled moments, rather than by a change of 
the exchange interaction or a destruction of the moments. 

C. Strong coupling: spin-charge interaction 

At large U, a quench within a Mott insulator freezes 
virtual charge fluctuations, leaving behind a certain den¬ 
sity ns of long-lived mobile carriers . The mechanism 
for the decay of the antiferromagnetic order is thus ex¬ 
pected to be the transfer of energy from excited quasi¬ 
particles to the spins, which is currently intensively in¬ 
vestigated in condensed matter pump-probe experiments. 
Although this mechanism is rather well understood in 
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FIG. 5: (Color online) Data for the following quench pro¬ 
tocol: t < 0 : J* = 0 (Neel state); 0 < t < 0.5: J* = 1, 
U = Ui\ t > 0.5: J* = 1, ?7 = 8. The intermediate step 
controls the excitation density in the final state, a) Time- 
evolution of the double occupancy or various values of Ui. b) 
Time-evolution of the order parameter M{t). c) Number of 
spin-flips per charge carrier density ns, [1 — M{t)]/ns, com¬ 
pared to the mean displacement R{t) of the initially localized 
particle in the t-J^ model for Jex = 0, 0.5 and Jex = 0 (dotted 
and dashed black lines, see text). The inset shows the the 
average of d{t) for 2 < t < 5 (open circles), and the density 
of mobile carriers ns, obtained from the integrated weight in 
the upper Hubbard band (red crosses). 


contrast to the dynamics at intermediate coupling, it is 
worthwhile to see how it can be investigated in the cold 
atom setup, because experiments in solids are very chal¬ 
lenging. 

To investigate the decay of long-range order systemat¬ 
ically, one has to vary the excitation density. Here we 
use a quench protocol where in addition to switching on 
the hopping at time t = 0, the interaction is changed 
to an intermediate interaction value Ui for a short time 
0 < t < 0.5, before it is set to the final value U for t > 0.5. 
(Note that various other protocols, such as an interme¬ 
diate time-dependent modulation of the hopping, would 
have the same effect.) Small values Ui lead to a larger 
double occupancy (Fig.[^), and indeed also a more rapid 
decay of M(t) (Fig.j^D). We also note that an exponential 


fit M{t) ^ -\- b would be consistent with a thresh¬ 

old behavior in which M (t) extrapolates to a finite value 
b for small excitation density {Ui close to Ui = 8) and to 
b — 0 for large excitation density, consistent with earlier 
quench studies based on the non-crossing approximation 
impurity solver [^ . but the times are not sufficient to 
analyze this long-time behavior in detail. 


For a quantitative analysis of the short time behav¬ 
ior, we determine the number ns of doublons and hole 
carriers in the quenched state (Fig. [St inset). Due to 
virtual charge fluctuations, ns is not exactly given by 
an instantaneous expectation value d{t) in the Hubbard 
model, and we compute ns from the total weight in the 
upper Hubbard band (Aon. IBl) [^ . For small times, the 
curves M{t) for various values Ui can then be scaled on 
top of each other by plotting [1 — M{t)]/ns (Fig. |St). 
Such a scaling implies that the number of flipped spins, 
1 — M{t), is proportional to the number of carriers. This 
is consistent with the picture that spin-flips are inserted 
by mobile carriers, which are initially localized and thus 
act independently up to times depending on ns- For large 
times there is a deviation from the scaling due to the 
gradual melting of the order parameter. 


To further corroborate this picture, we analytically 
compute the spin-flip rate per carrier in the low den¬ 
sity limit from the behavior of a single carrier which is 
initially localized at a given site 0 in a Ising spin back¬ 
ground. Following Ref. [2^, we omit the transverse dy¬ 
namics of the spins, which is on the timescale of Jex 
and much slower than the hopping, and keep only the 
2 -component of the exchange coupling Jex {t-Jz model). 
The model can then be reduced to a tight-binding model 
for a single particle on the lattice, with effective Hamil¬ 
tonian H = -J*l\fZY.{i 2 ) the 

number of flipped spins is simply given by the displace¬ 
ment |j| from the origin, and the second term in the 
Hamiltonian accounts for the corresponding exchange en¬ 
ergy cost, i.e., the particle is bound to the origin by a 
linear potential due to the “string” of flipped spins left 
behind [13 . The dotted line in Fig. (St shows the mean 
displacement R(t) of the particle in this model, which in¬ 
deed coincides with the mean number of flipped spins per 
particle in the numerical DMFT results. As evident from 
a comparison of the two curves for Jex = 0 and Jex = 0.5 
(the perturbative value for the Hubbard model at 17 = 8, 
see also Fig.|4^), the effect of Jex becomes important only 
at longer times (when numerical data already depend on 
ns), because initially the kinetic energy of the carrier is 
much larger than Jex. In order to measure the effect of 
Jex on the charge-carrier interaction, one would have to 
reduce the number of excitations (e.g., by switching on 
the hopping slowly), which however makes an accurate 
determination of ns increasingly difficult. 
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IV. CONCLUSION 

In conclusion, we have studied the short-time relax¬ 
ation dynamics of the Neel state in the single-band Hub¬ 
bard model by means of nonequilibrium DMFT, using 
DMRG to solve the quantum impurity model. We find 
qualitatively different relaxation behaviors for weak and 
strong interactions, separated by a crossover around U 
Ki 0.6X bandwidth: For strong interaction, local magnetic 
moments persist while their order is destroyed by spin- 
flips due to the hopping of mobile charges. The latter 
resembles the femtosecond carrier spin interaction which 
is relevant for th e dy namics of photo-induced states in 
high-Tc cuprates . To demonstrate the persistence of 
local moments we proposed a spin precession experiment, 
which could be implemented similar to the proposed mea¬ 
surement of dynamic spin-spin correlation functions in 
equilibrium (58| . At weak interaction, the dynamics of 
the Neel state is governed by almost conserved quasi¬ 
particles, which are also the origm for prethermalization 
in nearly integrable systems 0 , 1 , 0 - In the symmetry- 
broken state, the breakdown of these quasiparticles away 
from integrability leads to a crossover from oscillatory 
to non-oscillatory relaxation behavior, which can pro¬ 
vide a clear experimental signature that does now rely 
on a quantitative comparison to the thermal equilibrium 
state. 

Our simulations within DMFT are exact in the infinite 
dimensional limit, and it is thus interesting to compare 
to recent results for one dimension . Similar to our re¬ 
sults, in d = I one finds a rapid saturation of the double 
occupancy and a slower dynamics of the order parameter 
at large U, but the decay of antiferromagnetic order is of 
different origin: In large dimensions, the fastest melting 
processes after the quench take place on the timescale 
of the hopping due to the strong charge-spin interaction, 
while the latter is absent in d = 1 so that the dynam¬ 
ics happens on the timescale of the exchange interac¬ 
tion [^. The quasiparticle physics at weak coupling 
and in the crossover regime has not been addressed in 
Ref. [1^, but based on the perturbative argument given 
above the signatures in the off-diagonal components of 
the momentum distribution should persist also in lower 
dimensions. (Also in the paramagnetic case, a long-lived 
jump in the momentum distribution function is found in 

d =1 m, m, m and d =2 [n m.) 

Quench experiments starting from the Neel state have 
recently been performed with noninteracting fermions in 
one dimension [^ . and bosons in two-dimensions [2lj. 
Hence this setup should be a feasible approach to study 
fundamental aspects of the decay of antiferromagnetic 
long-range order in the paradigmatic Hubbard model. 
Moreover, on the numerical side our work emphasizes 
the high potential of DMRG as an impurity solver for 
future applications of nonequilibrium DMFT, to explore 
the intermediate coupling regime which is inaccessible by 
weak or strong coupling perturbation theory. 
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Appendix A: DMRG+DMFT setup 
General setup 

To simulate the dynamics of a lattice model which is 
initially in equilibrium at temperature T = 1//3, we adopt 
the formulation of dynamical mean-field theory within 
the Keldysh framework (nonequilibrium DMFT), for an 
L-shaped time contour C which extends from initial time 
t = 0 to a maximal time fmax along the real time axis, 
back to time 0, and along the imaginary time axis to —if3. 
For a general description of the formalism, as well as the 
notation and definition of contour-ordered functions, we 
refer to Ref. [d^. In this appendix we summarize the 
specific setup for the quench from the Neel state, and 
the solution of the DMFT equations using DMRG. 

In DMFT, the lattice model is mapped to a set of im¬ 
purity problems, one for each inequivalent lattice site j, 
with time-dependent hybridization functions 
The action of the impurity model is given by 

Sj “ y dtUnf{t)ni{t) 

-iY' [ dtidt2cl{ti)Aja{ti,t2)Ca{t2) (Al) 
a dc 

on the Keldysh contour C, which yields the lo¬ 
cal contour ordered Green’s function Gja{t,t') = 
Ca{t)cl{t')\/Z. The hybridization function 
Ajcr{t,t'), must be defined self-consistently. For the 
Bethe lattice, one has 

A,^{t,t') (A2) 

i 

where the sum runs over nearest neighbors of j. In the 
antiferromagnetic state, all sites on the A and B sub¬ 
lattices are equivalent, respectively. With the additional 
symmetry Ga,(t = Gb,-(t only one impurity model must 
be solved with A„{t,t') = where we 

used the scaling J{t) = with the coordination 

number Z. For the initial product state with J*(t) = 0 
for t < 0, A{t,t') = 0 if one time argument is on the 
imaginary branch of C. Furthermore equivalence under 
a simultaneous spin and particle-hole transformation im¬ 
plies the symmetry 

A>{t,t') = A<,{t,tr. (A3) 


To compute the Green’s function we follow Ref. 
and map the impurity model to a time-dependent An¬ 
derson Hamiltonian 

^^imp = U n-t-nj, -I- E ^pa0^p(jO>p(j +E [ fpcr (t flpcr + C . ] 

pa pa 

(A4) 

in which the impurity is coupled to L bath orbitals {p = 
The parameters Vpa-{t) and Cp are determined 
such that the local Green’s functions obtained from (EH) 
and (IA4I) are identical. As derived in Ref. [3 > for J* {t < 
0) = 0, one can choose Vpa{t) = 0 for t < 0, and the 
mapping condition is satisfied by (assuming L even) 

L/2 

A<{t,t') = zY,Vpa{t)v;At') (A5) 

p^l 

L 

A>(t,0 = -* E VpAtWplit'), (A6) 

p=L/2+l 

where Cpa- = 0, and the bath orbitals p = 1...L/2 and 
p = L/2 + 1...L are initially doubly occupied and empty, 
respectively. Equations (IA5|) and (IA6I) are solved by 
a Gholesky fit of the real-time matrix which 

quickly converges for small times with the number of bath 
orbitals required [i^. Due to the symmetry (IA3I) we use 

Vp^.^it) = VL/2+p,a{t)* for P < L/2. (A7) 

The impurity site is initially occupied with a spin a =t 
(for a site on the A sublattice), i.e., the initial state 
for the impurity model is a product state I'kimp^^) = 

c]-rifci Green’s function is obtained 

by solving 

t') = *(^imp,A|ct(i')c^ Wl^imp.A), (AS) 

GXai^A') = -i(^'imp.A|Ca(t)cj,(t')l^imp.A), (A9) 

where time evolution is determined by (lAill . We use a 
Krylov time-propagation for matrix product states 
with up to L = 24 bath orbitals. 


Inhomogeneous setup 

For the inhomogeneous setup we assume that in the 
initial state on the lattice the spin at one site o of the 
lattice is flipped in the x direction. Without loss of gen¬ 
erality we assume that o is on the A sub-lattice. From 
the self-consistency equation (IA2I) one can see that the 
hybridization on all other sites differs from the homo¬ 
geneous case only in order 1/.Z, i.e. for Z —^ oo the 
back-action of the probe site on the rest of the lattice 
can be neglected. On the probe site we solve an impurity 
problem with the same (nonequilibrium) hybridization 
function Ay4 as on all remaining A-sites, i.e., an impurity 


problem (IA4I) wit the same parameters Vpo-, but with a 
different initial state, 

Lj2 L/2 

I'l'imp.o) = apiUpj^lO) -f-cj ap.|apj^|0)^/-\/2. 

i=l i=l 

(AlO) 

Observables 

Local observables = (4>inipj|(!I(t)|'kimpj) are 

directly measured in the impurity model (j = o. A), in 
particular the density O = Ua-, the double occupancy 
O = rifni, and the spin O = Sa=x,y,z = ^ J2aa' 4'''aCo-' 
(tq, are the Pauli matrices). 

In the translationally invariant case (no probe site), we 
also determine diagonal and off-diagonal components of 
the momentum occupations n(e,t) and 7r(e,t), which are 
obtained from the momentum resolved Green’s function 
(for the definition of k and fc, see the main text) 

r (f Tl = E*^^ccfe(t)cj;,(<')) -i{Tcck{t)cl{t'))\ 
\-i{TcCk(t)clit')) -i{TcCi,it)cl(t')}j ' 

(All) 

(Here and in the following, bold-face quantities denote 
2x2 matrices and we omit spin indices for simplicity). 
The self-energy is local in space but depends on the sub¬ 
lattice and spin; in the k,k representation it thus assumes 
the (2 X 2) form 

S(t,t') = i[EA(t,t') + SB(t,t')]l 

+ i[EA(<,t')-SB(t,t')]^x, (A12) 

so that Ge is obtained from the lattice Dyson equation 
Ge = {idt + p — e — 'E)~^, where the dispersion in the 
k,k representation reads e = er^ because = —e^. The 
components T,j of the self-energy (j = A, B) are obtained 
from the impurity Dyson equation {idt+p—Aj — 'Etj)~^ = 
Gj. In praxis, we solve an integral equation Gj = Zj + 
Zj * Aj * Gj for Zj = (idt + p — Sj)~^. We then have Z = 
(idt+p—Tj) ^ = ^(Zji + Zs)l-h^(Zy^ — Zs)T^, and G^ is 
obtained from the integral equation G^ = Z + Z * e* G^- 

Appendix B: Mobile carrier density in the excited 
state. 

In the Mott insulating phase of the Hubbard model, 
a well-defined measure for the number of doublon or 
hole carriers is given by the total occupied spectral 
weight in the upper Hubbard band and the total unoc¬ 
cupied weight in the lower Hubbard band, respectively. 
The double occupancy, in contrast, depends on virtual 
charge fluctuations which are nonzero also in the insu¬ 
lating ground state. Specifically, we define the occu¬ 
pied density of states as the partial Fourier transform 
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FIG. 6: a) Occupied density of states Na{uj,t) after the 
quench t < 0 : u = 0 (Neel state); 0 <t < 0.5: v = 1, U = Ui; 
t > 0.5: u = 1, [/ = 8; for various values of Ui] and 

N_i{oj,t) refer to majority and minority spin, respectively, and 
the upper Hubbard band is scaled by a factor 5. Dashed and 
solid lines denote time t = 4 and the largest simulation time 
t = tmax, respectively, b) Integrated weight in the upper 
Hubbard band (red crosses and blue stars). 


chain by introducing operators which are invariant un¬ 
der all permutations of the branches of the Bethe lattice 


Cn 




(Cl) 


where Zn = X)r|i-o|=n number of sites on the n-th 

nearest neighbor shell. Then the action of the Hamilto¬ 
nian is determined by [H, Cj] = — with 


h = 


/O 1 0 0 --A 
10 10 
0 10 1 




(C2) 


Hence, eigenvectors for the eigenvalue e satisfy the equa¬ 
tion 


Nc,{t,uj) = Im Jp cisexp(—s^/2(5^) exp[—isa;]G^(t — s,t) 
where G^{s, s') = i{c\{s')c„{s)) is the local Green’s func¬ 
tion, and 5 = 1.5 ensures a smooth cutoff (which does not 
influence the results unless its inverse width is longer than 
the inverse of the gap). The spectrum N„{uj, t) is plotted 
in Fig. [6^ for two different times, for the same quench 
parameters as in Fig. Oof the main text. The right panel 
shows the integrated density WA(t) = dujN„{u},t). 

While the weight in the upper and lower band differs 
considerably between majority and minority spin, the in¬ 
tegrated weight Wa{t) reflects the doublon density and 
is thus independent of a. It is interesting to point out 
that as a function of time spectral weight is both redis¬ 
tributed between the lower Hubbard bands of the two 
spin components (which reflects the decay of the Neel 
order), and within the upper Hubbard band (which re¬ 
flects the change of the kinetic energy of the doublons), 
while the total weight in the upper band is roughly con¬ 
stant (see, e.g., Ref. [3). For the analysis in the main 
text, we take ns = (WA(^max) + WA(imax))/2. 


A(e)o = 1 

(C3) 

A(e)i = e 

(C4) 

eA(e)n = A(e)n+i + A(e)n-i, 

(C5) 


and are thus given by the Chebychev polynomials of sec¬ 
ond kind [^, 4'{e)n = [/„(e/2) for —1 < e/2 < 1. The 
Un can be conveniently written as 


Un{cos{9j) 


sin[(n -I- 1)9] 
sin(0) 


(C6) 


from which one can also see the orthogonality 


dxw{x) U„{x)Um{x) = Sr, 


(C7) 


with w{x) = ^\/l — x^. Thus the solution of the Heisen¬ 
berg equations of motion for the local c operator m 


Appendix C: Solution for U — 0 

For U = 0 the time evolution of the Neel state on 
the Bethe lattice can be obtained analytically by solving 
the Heisenberg equations of motion for the c-operators, 
which provides a good check for the numerical implemen¬ 
tation. For completeness, we provide this solution in the 
following. We choose site 0 to be the origin of the Bethe 
lattice, which is on the A sub-lattice without loss of gen¬ 
erality. One can map the solution of equations of motion 
on the Bethe lattice to a one-dimensional semi-infinite 


^Go(t)=f[H,Go(t)], Go(0)=Go, (C8) 

is given by 

OO 

Co{t) = J2Mt)Cr„ (C9) 

n—0 

An = y dxw{x) e~'‘^'''*Un{x). (CIO) 
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This can be transformed to 


= - /' dicos{9)) sin(0)e-*2 

TT sin^t*) 

= - / d6»sin[(n + 

TT _/q 2 zt 




d6»e-*2co«(»)‘agsin[(n + 1)6»] 


z(n + 1) 


irt 


cos[{n + 1)9] 


!h±3. f d6»e*2™'^(«)*cos[(n +l)6l](-l)"+i 
7 ^^ Jo 


= Hr(n + l) 


>/n+l( 2 t) 


(Cll) 


The second to last line is a variable transformation 
0 —>■ 71 — 9, and in the last line we have used the integral 
representation of the Bessel function M, 

= d0cos(n6l)e*^“"(®). (C12) 

TT Jo 

The explicit form of the C-operators can be used to 
obtain local observables 

(Vf-lct (t)c.Ki)l^') =J2rnit)J^m{t){'f\ClCm.'m 

n,m 

(C13) 

= (C14) 

n 

where the expectation values are simple initial state val¬ 
ues. We start by evaluating the time evolution of the 
magnetic order, — nj,, at site 0, in the classical Neel 
state. For the latter we have 

(4'Neel|Cl.^C„t - = (-1)", (C15) 

and hence 

(^fNeellnf-it) - ni(t)\^mei} = 2^(-l) - 


n—0 




(C16) 


(where the summation index has been shifted by one). 
We can now use Gegenbauers addition theorem for Bessel 
functions [blj to obtain the final result 


('I'Neel|%(i) — 'ra4,(^)l'I^Neel) 


Jl(4t) 

2t 


(C17) 


which fits the numerics. 

Next we compute site 0 expectation values on the probe 
site. Now the initial state is a superposition 


\^) = (l^-Neel, t) + l^-Neel, i)) / ^2 (CIS) 

I^NeehO") = cJ_^Co,-f-|4'Neel)- (C19) 

We evaluate the cross-spin expectation values 

('5'c}'(t)) = {'l>\cl^(t)coi{t)'5>) 

= (C20) 

n 


Spin-flip expectation values are only non-zero in the ini¬ 
tial state at site 0, where we have 


{S^t)) = |V'o(OP(^l4coil'I') 

^ Ji(2t)" 

2^2 ■ 


(C21) 

(C22) 


Hence S';)" (t) is purely real, so that the dynamics is en¬ 
tirely longitudinal in the Sx-Sy-plane, 


(^0 W) = (C23) 

{Sm = 0- (C24) 
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